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Abstract 

The complex magnetic field of the 
permanent-magnet electron cyclotron 
resonance (ECR) ion source MEFISTO 
located at the University of Berne have 
been numerically simulated. For the first 
time the magnetized volume qualified for 
electron cyclotron resonance at 2.45 GHz 
and 87.5 mT has been analyzed in highly 
detailed 3D simulations with unprecedented 
resolution. New results were obtained from 
the numerical simulation of 25211 electron 
trajectories. The evident characteristic 
ion sputtering trident of hexapole confined 
ECR sources has been identified with the 
field and electron trajectory distribution. 
Furthermore, unexpected long electron 
trajectory lifetimes were found. 

1 Introduction 

Unlike the magnetic field of solenoid coils the field 
of permanent magnets can only be tuned over a 
limited range after the manufacturing process and 
the construction of the ECR ion source. Great care 
has therefore to be taken in the design of the mag- 
netic configuration because empirical adaptation of 
the magnetic field after construction is limited to 
very few options. The magnetic arrangement cho- 
sen for the MEFISTO [4 ion source features a field 
distribution which cannot be reduced from three to 
two dimensions for the purpose of a more economic 
simulation. Hence appropriate 3D simulations have 
to be performed to bring out the full spatial extent 



of the field distribution. 




Figure 1: The arrangement of the full permanent 
magnets of the MEFISTO ECR ion source at the 
University of Bern (Marti et al, 2001). Some mag- 
nets are removed in this drawing for better visual- 
ization 

Figure [I] shows the arrangement of the perma- 
nent magnets in the MEFISTO ion source. The six 
permanent magnets of the ring-shaped cluster on 
the left are all magnetized inwards. Similarly the 
cluster on the right side is magnetized outwards. 
These two outer ring-shaped clusters form the mag- 
netic field of an axial magnetic bottle and thereby 
establish the axial confinement of the plasma elec- 
trons. In contrast the central magnet cluster con- 
sists of permanent magnets which are magnetized 
alternatively inwards and outwards. This central 
cluster forms the hexapole field which produces the 
radial confinement of the plasma electrons. The 
ECR plasma is located inside this central ring ar- 
rangement. In addition Figure [I] also shows the 
helix antenna which radiates the microwave power 
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into the ECR plasma. 

The ionization process inside an ECR ion source 
works by electron collision with neutrals and ions. 
The electrons are effectively heated by the mi- 
crowaves, most effectively for those meeting the 
electron cyclotron resonance condition. Equation 
[I] shows that the ECR condition in general and the 
resonance frequency uo in particular depend on the 
local magnetic field density B, the electron mass 
m e and the charge of the electron e. However, the 
resonance frequency uj does not depend on the elec- 
tron energy as long as the electron is not acceler- 
ated close to the speed of light where it's change of 
mass would come into account due to special rela- 
tivity. 



^ECR 



eB 



(1) 



Therefore the choice of the excitation frequency 
band determines the required magnetic field den- 
sity B. Hence the design parameters of the magnet 
structure have to be chosen carefully so that the 
magnetic field fulfills the ECR criteria 0. This cri- 
teria is generally fulfilled on a closed iso surface of 
the magnetic field density (see Section [3]). This iso 
surface defined by the magnetic field distribution 
and the microwave frequency has to be located en- 
tirely inside the plasma container for optimal ECR 
efficiency. 

This ECR ion source is used for the calibration 
of space craft instruments which are designed to 
detect highly charged ions in space. The perfor- 
mance of the ECR ion source with respect to the 
production of these highly charged ions depends 
on the square of the ECR resonance frequency [2]. 
In addition the magnetic field distribution forms a 
magnetic mirror with a ratio r m = B/B centre in all 
directions, which should be as high as possible for 
optimal electron and ion confinement. However the 
successful extraction of highly charged ions requires 
a locally poor confinement performance so that the 
ions can be extracted from the loss cone in the axial 
beam direction. In Figure [I] this extraction system 
is formed by the two cylinder shapes to the right. 

All of these requirements have to be met for 
the successful operation of an ECR ion source. 
To achieve these goals three dimensional numerical 
simulations of the magnetic field can be used which 
allow the highly accurate design of future ECR ion 



sources. The simulation technology has been tested 
and the results have been compared with measure- 
ments of the existing MEFISTO ECR ion source. 

Due to the unique field distribution of hexapole 
confined ECR ion sources the high energy electron 
population features a trident shaped spatial dis- 
tribution as predicted theoretically by numerous 
authors (PQ, [2]) and verified in experiments by 
Zschornack et al. ([3]). The presented magnetic 
field simulation confirms both theory and experi- 
ment in three dimensions. 

2 The finite element model 

To achieve the highest possible simulation preci- 
sion given the available infrastructure a finite el- 
ement solver has been chosen. Magnetic simula- 
tions have been performed using Tosca by Vector 
Fields (www.vectorfields.com). All pieces and vol- 
umes relevant for the simulation of the stationary 
magnetic field were modeled carefully. The result- 
ing finite element mesh itself has been tested exten- 
sively with reference patterns and circular charged 
particle trajectories such as electrons and ions. No 
deviation from the circular orbits were discovered 
within the system limit of 5000 integration steps 
for each particle. A summary of the FEM mesh 
characteristics is given in Table [I] 



number of linear elements 
number of quadratic elements 
number of nodes 
number of edges 
RMS 

distortion of worst element 



2.4* 10 6 
1.9* 10 5 
7.2* 10 5 
3* 10 5 
1.033 
8.31 * 10- 4 



Table 1: Some characteristics of the finite element 
model. 

Figure [2] shows the finite element model used to 
simulate the magnetic field. 

As described in Section [T] the model consists 
of three sets of block magnets arranged in rings 
around the central plasma area and two massive 
soft iron rings separating each of the block magnet 
ring arrangements. Figure [3] shows the magnetic 
field density along the beam axis with respect to 
the z-coordinate. The ECR zone is indicated at 
87.5 mT where the field density meets the ECR 
condition. 
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surement. All subsequent simulations and results 
are based on the results of the presented numerical 
field simulation. 



3 The ECR zone 

The ECR ion source is fed with a constant mi- 
crowave frequency of 2.45 GHz, so the location of 
the resonance region can be determined by the cor- 
responding magnetic field density B res = 87.5 mT 
0. To present a two dimensional analysis of the 
field distribution and the ECR locality inside the 
magnet system we introduce cut planes (Plane A, 
B, C) into the three dimensional model as shown 
in Figure [4] 



Figure 2: The finite element model with its essential 
parts. 




-200 -150 -100 —50 50 100 150 200 



Position [mm] 




Figure 4: To map the magnetic field density onto 
two dimensional patches we introduce cut planes in 
the 3D model. 



Figure 3: Simulated magnetic field density along 
the beam axis. 



The axial magnetic field of the existing 
MEFISTO ECR ion source has been measured pre- 
viously [4] [6]. The local maxima and the local 
minimum of the measured magnetic field is given 
with 240 mT and 60 mT respectively. The corre- 
sponding mirror ratio is r m = 240 mT/60 mT = 4. 
The maxima of the numerical simulation are 240.1 
mT, the minimum is 60.05 mT and the mirror ratio 
r m = B max / B m i n is 3.998. The numerical model 
is therefore in excellent agreement with the mea- 



In the following figures the magnetic field density 
has been mapped to color space. There are two 
figures of each cut plane (see Fig. [5] to 10) giving 
an overview and a detailed view of the central part 
of each plane respectively. Iso contour lines of the 
field density are also indicated on each graph. 

Because cut plane B is rotated by 90 degrees with 
respect to cut plane A the magnetic field is cut at 
a different angle too revealing its non symmetric 
distribution in cut plane B (Figure [8]). This feature 
will be visible more clearly in the 3D view presented 
further down this Section. 

Despite the hexapole arrangement of the central 
magnet cluster the resulting field distribution in the 
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Figure 5: Simulated magnetic field density in plane 
A. 




Figure 6: Detailed view of the center part of plane 
A. The ECR iso contour line is highlighted. Due 
to the orientation of the cut plane (see Fig. ^ the 
shape of the iso contour line is almost symmetric. 

symmetry plane C (see Fig. [4| yealds concentric iso 
contour lines as presented in Figures [lOj However 
this is true only for contour lines which are close to 
the center of the hexapole. Contour lines further 
from the center and closer to the actual block mag- 
nets show the expected hexapole pattern with its 
distinct poles. 

The small irregularities within the contour lines 
are due to the discretization of the model by the 
finite element size. The resulting relative error in 
the field density is lower than 10 ~ 2 and has been 
chosen to be acceptable. 
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Figure 7: Simulated magnetic field density in plane 
B. 




Figure 8: Detailed view of the center part of plane 
B. The ECR iso contour line is highlighted. Cut 
plane B is rotated by 90 degrees with respect to 
plane A. The contour lines are no longer symmetric 
with respect to the beam axis. 

In each detailed view above the iso contour line of 
87.5mT is highlighted. Extending the iso contour 
line concept from two to three dimensions leads to 
an iso contour surface. This iso contour surface 
includes all element centroids of the finite element 
model which fulfill the criteria of the given mag- 
netic field density B re f = 87.5 mT. On this sur- 
face the gyration frequency of the electrons in the 
magnetic field equals the microwave frequency. The 
plasma electrons are therefore heated effectively on 
this surface. The surface defined by this criteria is 
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Figure 9: Simulated magnetic field density in plane 

c. 
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Figure 10: Detailed view of the center part of plane 
C. The ECR iso contour line is highlighted. To- 
ward the center the iso contour lines are concentric 
circles. 



shown in Figure 11 



On Figure [TTJ the centroid of each FEM element 
of the specified surface is represented by a small 
bubble. Each bubble also represents one launch 
coordinate for the trajectory model represented in 
the following section. The shape is symmetric with 
respect to the origin located at the center of the 
shape. Note the trident shapes on both ends. These 
tridents can be explained as follows. Because the 
locations of the trident shapes have an offset from 
the center of the magnet arrangement in positive 
and in negative z direction the radial magnetic field 




Figure 11: Isosurface of a constant magnetic field 
density B res (upper image) and cut-away 



component of the outer ring magnet clusters no 
longer cancel out. This leads to a dominance of 
the outwards magnetized cluster in positive z direc- 
tion. Consequently the radial field of the outward 
magnetized hexapole magnets is enhanced and the 
field of the inwards magnetized hexapole magnets 
is suppressed where the trident shape in positive 
z direction is located. The trident shape in neg- 
ative z direction is formed in the same way with 
the field of the inward magnetized hexapole mag- 
nets suppressed and the field of the outwards mag- 
netized hexapole magnets enhanced. Hence, the 
shapes are twisted relative to each other by 60 de- 
grees because of the same angle separating inward 
and outward magnetized permanent magnets in the 
hexapole cluster. The cross section through the 
mid plane is circular (see Fig. 



10 for 2D contour 



plot) and has a diameter of 21.0 mm. The overall 
length of the shape is 37.3 mm. 

The described trident shapes of the ECR zone 
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also influence the plasma as explained further 



down. Figure [12] shows a photograph of two super- 
imposed sputtering tridents observed in a hexapole 
confined ECR facility. 
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Figure 12: The ion beam of hexapole confined ECR 
ion sources produces trident shaped sputtering pat- 
terns (Picture from CASYMS ECR ion source, 
University of Bern). Here two trident patterns are 
superimposed on a microwave antenna. This an- 
tenna was used to radiate the microwaves into the 
plasma and for impedance matching. It was in con- 
tact with the ECR zone. 

Cut plane C is located at z=0 which results in cir- 
cular iso contour lines of the field density near the 
center not similar to the observed trident shaped 
sputtering patterns. However if the cut plane is 
shifted parallel to z=17mm the iso contour lines 
be identified with the sputtering trident (see Fig. 
13). This is because the shifted contour plane cuts 



the isosurface (Fig. 11) no longer at the center (as 
plane C is doing it) but 17mm toward the trident 
face. 

We can clearly identify the top face of the shape 
with the typical sputtering trident well known for 
ECR ion sources with a hexapole radial confine- 
ment such as MEFISTO. The identification of the 
simulated and the observed trident takes into ac- 
count the shape of the ECR region, the presence of 
high energetic electrons and the subsequent emer- 
gence of ions as follows. 

According to the successful ECR plasma model 
of Wurz et al. ([12 j) the ionization is a step by 
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Figure 13: Simulated magnetic field density in a 
plane parallel to C, but shifted by 17mm. While 
contour plot C at z=0 (Fig. [7fl[ ) features circular 
iso contour lines another contour plot at z=17mm 
reveals the similarity to the sputtering trident. The 
three separated islands of the ECR iso contour lines 
are highlighted. 



step process starting with few free electrons and a 
majority of neutrals. The few free electrons get ac- 
celerated to high energies by the electric field of the 
microwave injection in the described way of electron 
cyclotron resonance. The energetic electrons then 
collide with neutrals to produce single charged ions 
and more free electrons. More neutrals get ionized 
due to the higher electron density. Simultaneously 
single charged ions collide with energetic electrons 
to be ionized from single to double charged ions and 
so forth. This way very few free electrons create an 
avalanche effect and finally establish the quasi sta- 
ble electron density of the ECR plasma featuring 
highly charged ions. A fully consistent version of 
this model has been tested numerically and verified 
up to a charge state of Ar8+ by Hohl et al. ([13]). 

The magnetic trident face has been identified as 
part of the volume in which ions are created. Ions 
are also created by the same process in the rest of 



the depicted shape (Fig. 11) as well as inside of 
it. However only the trident part of the front face 
is being mapped by the extraction ion optics to 
the target surfaces where the sputtering has been 
observed. 

The magnetic field density increases in every 
direction pointing away from the entire depicted 
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shape. The closed depicted shape therefore con- 
fines a volume of minimal magnetic flux density. 
The magnetic arrangement thereby fulfills the cri- 
teria of a minimum B-field structure. Further we 
have the liberty to refer to the defined volume as 
a plasmoid - a plasma magnetic entity [9]. The 
presented isosurface element centroids were used to 
define the launch points of the electron trajectories 
described in the next section. 



4 The trajectory model 

The ECR plasma is charged slightly negative as 
described by Shirkov [10], [11] and verified experi- 
mentally by Golovanivsky and Melin [7 j. Accord- 
ing to these measurements the plasma potential af- 
fecting the electron trajectories has been chosen to 
be minus 10 V. Due to the compact shape of the 
simulated plasmoid the spatial distribution of the 
plasma potential was modeled as a sphere with the 
same diameter as the plasmoid. 

The electron heating takes place within the ECR 
qualified volume sheath with a magnetic field den- 
sity of 87.5 mT including an assumed tolerance of 
±0.5 mT, which is due to the bandwidth of the 
implemented microwave generator. Because of col- 
lisions and momentum transfer the production of 
new ions can only happen close to where the elec- 
trons are energized. Consequently the initial launch 
coordinates of the simulated electrons were cho- 
sen to be located within the ECR qualified volume 
sheath. Each bubble seen in Figure [TT] is located 
at an element centroid of the magnetic FEM model 
and therefore represents one launch point for the 
trajectory model. In total the iso surface shape 
yields 25211 launch points. 

The ECR heating process mainly takes place 
perpendicular to the local magnetic field lines [2]. 
Therefore the initial velocity distribution has been 
chosen to be anisotropic. The velocity component 
parallel to the local magnetic field of each trajec- 
tory was scattered symmetrically around zero with 
a maxwellian temperature of 2eV [8]. Figure 14 
shows a histogram of the modulus of this vector 
component. 

In contrast the modulus of the two velocity com- 
ponents perpendicular to the magnetic field were 
scattered around a kinetic particle energy of 2keV 
corresponding to the measurements of M. Hohl et 




1.0 1.5 2v0 
Velocity^]' 



0.U 



Figure 14: The velocity component parallel to the 
local magnetic field is distributed with a correspond- 
ing temperature of 2e V. 



al. [5]. Figure 15 shows a histogram of the modulus 



of the perpendicular velocity components centered 
around 2keV with a maxwellian temperature of 2eV 
with respect to the offset. 
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Figure 15: The norm of the velocity component per- 
pendicular to the magnetic field has an offset of 
2ke V and a corresponding temperature of 2e V with 
respect to the offset. 

The model does not take into account any par- 
ticle scattering phenomena nor any other kind of 
particle particle interaction. The element size in 
the finite element and the selection criteria of the 
ECR zone gave a total number of launched parti- 
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cles of 25211. All simulations have been performed 
on a IntelCore2Duo 6700 clocked at 2.66GHz and 
2 GB of RAM. The average run time was 83 hours 
at more than 99 % of CPU load. 



simulation will be referred to as simulation A. Fig- 
ure 17 gives the mean velocity of each trajectory 



with respect to its length for this simulation. 
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5 Results 

The trajectory lifetimes, lengths and velocities of 
each particle were analyzed. Due to hardware limi- 
tations the simulations were stopped once less than 
one percent of the initial number of electrons re- 
mained in the simulation. The life time of a par- 
ticle is limited by the collision with a wall of the 
plasma container or by the end of the simulation. 
In a histogram of electron trajectory counts versus 
trajectory life time we can fit two exponential decay 
functions of the electron life times t\ and T2: 



t.\ ( --Mn(l/2) . 1 --Mn(l/2))\ / x 

n(t) = no ( e T i w ; + —e T 2 w 1 (2) 

Figure [16] shows a histogram of the particle life 
time and an exponential decay with r± = 10/is and 
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Figure 16: Histogram of particle life time and the 
corresponding exponential fits. 

There are two significant deviations from the ex- 
ponential decay behavior around 136/is and lesser 
so at 242 11s. The deviation represents two electron 
fractions featuring unexpected long life times. This 
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Figure 17: Mean velocity versus trajectory length of 
simulation A. Two clusters at 3537m and 6287m 
are separated from the main population 

Apart from the main population from to 2000m 
there are two clusters of trajectories accumulating 
at a trajectory length of 3537m and another one 
at 6287m. Both cluster areas present a wide ve- 
locity spectrum centered around the same value as 
the main population. The corresponding trajecto- 
ries have been traced back to find their underlying 
launch points. The respective launch points and 
their initial conditions are distributed both homo- 
geneous over the plasmoid shape as well as correctly 
anisotropic in velocity space. Hence, the phenom- 
ena is not the result of a special launch location of 
the trajectories nor is it due to an anomaly of the 
initial velocities. 

To distinguish the simulated phenomena from 
numerical errors associated with the random par- 
ticle velocity distribution a second simulation was 
performed. The second simulation made use of the 
same principle of anisotropic initial velocity distri- 
bution as the first simulation but starts with a new 
set of particles. Due to the randomized velocity 
vector creation process the new set is deliberately 
different from the first one despite the equal mean 
energy and the equal initial start locations. Figure 
18 shows the mean velocity versus trajectory length 
of the second simulation B. 
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Figure 18: Trajectory length versus mean velocity 
of electrons in second simulation run B. 

Table [2] gives a statistical breakdown of the tra- 
jectory lengths of both simulations. 

Simulation \i a 

First simulation A 296.1m 1207.6m 

Second simulation B 343.8m 1213.5m 

Table 2: A statistical breakdown of the trajectory 
length of both simulations, n refers to the median 
and a to the standard deviation. 

The high standard deviation indicates the wide 
spread of the length distribution. The medians of 
both simulations show a shift of 47.7m. We can see 
the graphs for both simulation runs are similar in 
overall shape. Both show the same two clusters in 
the same area. However in the second simulation 
the clusters are slightly shifted toward a higher tra- 
jectory length which corresponds to the higher me- 
dian of the second simulation. Table [3] summarizes 
the findings related to both clusters. 

Due to the different initial conditions of velocity 
distribution the two simulation runs are not iden- 
tical. Both clusters show a relative shift between 
2.59% and 3.26%. However both simulation runs 
show the clusters are clearly separated from the 
main population and grouped in the same range of 
trajectory length. Neither median of either sim- 
ulation corresponds to a dimension given by the 
plasma container nor a low multiple integer of such. 



Simulation l^clusterl I^cluster2 

First simulation A 3537m 6287m 
Second simulation B 3631m 6499m 
A rel 2.59% 3.26% 

Table 3: A summary of the clusters separated from 
the main population. Hdusteri refers to as the me- 
dian of the first cluster in each simulation and 
I^duster2 to the second one. Arel depicts the rel- 
ative difference between each median. 

6 Discussion 

Moving electrons in a magnetic field are deviated 
from straight trajectories by the Lorentz force per- 
pendicular to their flight directions. The electrons 
enter circular orbits with a radius r^. This is called 
the Larmour radius of an electron. It depends 
solely on the electrons energy and the local mag- 
netic field. It is given by: 



where v is the velocity of the electron, m its mass, 
e its charge and B is the local magnetic field den- 
sity. For the case of the MEFISTO ECR zone and 
2keV electrons this results in = 1.722mm. The 
Larmor radii in the ECR zone are therefore much 
smaller than the plasma container. Hence we can 
give a rough approximation of revolutions an elec- 
trons takes before colliding with a container wall 
with fi as the median of the first simulation run 




n rev ^ /i/(27rr L ) = 2.7 * 10 4 (4) 

Given a collisionless plasma the electron motion 
in the simulated magnetic arrangement is there- 
fore dominated by the magnetic field rather than 
the container dimensions and could indeed produce 
long trajectories as it is suggested by the simulation 
results. However this does not apply for the current 
operation parameters of the MEFISTO ion source 
as we will see in Section \6J\ 

6.1 Relevance to laboratory plasmas 

Laboratory plasmas as in the MEFISTO facility are 
not collisionless. The mean free path is limited not 
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only by the plasma container walls but also by colli- 
sions with other particles of the plasma. According 
to Chen [1 the mean free path of an electron is 
given by 



A = 3.4* 10 



li 



nln A 



m 



A 127mA!) 



Ad is the Debye length given by: 



Ai 



e KT e 



(5) 



(6) 



(7) 



where T e y is the electron temperature, K is the 
Boltzmann constant, n is the plasma density and e 
is the charge of the electron. Assuming a plasma 
density of 10 16 m -3 and an electron temperature 
of 2keV we get a Debye length A^ = 3.24mm, a 
A of 1.278 * 10 10 and a mean free path length A of 
5.84m. 

This value clearly shows the mean free path 
length of 2keV electrons in an ECR laboratory 
plasma are far too short to produce the phenom- 
ena suggested by the trajectory simulation. How- 
ever a plasma density of 10 13 m -3 would lead to 
a mean free path of 5089m by the same calcula- 
tion. This would bring the simulation results into 
physical possibility. 

If we assume an effective cross section for electron 
collisions with neutrals of 10 -20 m 2 [2] and again 



a plasma density of 10 



16 



m 



we get a mean free 



path to the next neutral particle of A n = 10 4 m. 
Collisions with neutrals are therefore not limiting. 

6.2 Relativistic considerations 

Due to the electron energies involved of 2keV a 
slight change of electron mass occurs due to special 
relativity. The electron cyclotron resonance condi- 
tion is attached to the electron mass and the mag- 
netic field density 0. A change in electron mass 
therefore results in a different resonance condition 
of the qualified magnetic field density. The relation 
is as follows: 



Bi 
Bo 



mi 
m 



1 



Where v is the velocity of the electron, c the 
speed of light in vacuum, B\ the ECR qualified 
magnetic field density with the consideration of 
mass increase and Bo without the consideration of 
mass increase, v/c equals 8.67* 10 -2 and Bi/Bq — 1 
equals 3.78* 10~ 3 . This modified magnetic field re- 
quirement would be fulfilled at a different location 
due to the spatial distribution of the magnetic field. 
In fact the isosurface depicted in Section 3 would be 
enlarged toward the plasma chamber walls where 
the magnetic field features a higher modulus and 
where the resonance condition would again be ful- 
filled. However the shift of the isosurface would be 
smaller than the thickness of the surface due to the 
implemented frequency tolerance of the microwave 
generator. Hence this drift in location of the elec- 
tron cyclotron resonance toward a higher magnetic 
field locality due to the effect of special relativity 
has been neglected. 



7 Conclusions 

The plasmoid shape resulting from the magnetic 
configuration of MEFISTO is responsible for the 
typical sputtering geometry of hexapole confined 
ECR ion sources. Future permanent magnet ECR 
ion sources can be designed very accurately using 
finite element modeling of the magnetic field. The 
accuracy of the FEM model was confirmed by com- 
paring it to experimental data on MEFISTO. The 
location of the ECR qualified plasmoid constrained 
by an isosurface of constant magnetic field density 
depends on the magnetic arrangement only. It can 
therefore be chosen arbitrarily by tuning the design 
parameters of the magnetic setup. 

Confined by permanent magnet hexapoles, long 
electron life times can be expected for low density 
plasmas in ECR ion source. This could also be very 
useful for electron or positron traps or any other 
particle storage facility as the numerical model is 
scalable by particle mass, particle energy and mag- 
netic field density. 
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